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When electron correlations are important it is often neces- 
sary to use numerical methods to solve the Hamiltonian for a 
finite system (cluster) "exactly". Unfortunately, such meth- 
ods are restricted to small systems. We propose to combine 
the "exact" numerical diagonalization for small clusters with 
the perturbative calculations to take into account the intra- 
cluster as well as intercluster interactions. 



I. INTRODUCTION 

Since the discovery of high-T c superconductivity^, the 
behavior of strongly correlated electronic systems re- 
mains a central problem in contemporary condensed mat- 
ter physics. In spite of considerable effort devoted to the 
analysis of these systems, it is clear that the necessary 
theoretical skills and tools to deal with strongly corre- 
lated fermion systems are lacking. There are no exact 
solutions except in one dimension (e.g., the t — J modelis 
exactly solved by the Bethe-ansatz method for J = 2£o) 
and approximate analytic techniques have been known 
to lead to qualitatively incorrect predictions. The fun- 
damental obstacle which appears in the analytical ap- 
proaches is the difficulty in handling the strong corre- 
lations in a satisfactory way. Moreover, in mean-field 
calculations it is necessary to make a priori assumptions 
about the ground-state properties. 

Therefore, most work on models with strongly corre- 
lated electrons has been done using numerical techniques. 
Among others, variational calculations!], various realiza- 
tion of quantum Monte Carlo simulations p and an exact 
diagonalization of small systems! are used to obtain the 
properties of these models. 

Unfortunately, numerical methods also meet some seri- 
ous problems. The main difficulty in the quantum Monte 
Carlo calculations is the famous sign problem, which re- 
duces the usage of this method at low temperatures and 
at the physically interesting densities. The minus sign 
problem does not arise in the diagonalization procedures, 
based on the Lanczos methodEl and its modifications^, 
where all quantities (static and dynamical) can be com- 
puted from the ground state. Regretfully, the Lanczos 
technique is limited to small clusters by the rapid in- 
crease of the size of the Hilbert space with the number 
sites. Typically the calculations are performed on 4 x 4 
cluster with periodic boundary conditions for one hole, 
two holes, or an arbitrary number of holes. With respect 



to the infinite lattice this corresponds to an investigation 
of only a small number of points in the Brillouin zone. 
Therefore, the overall shapes of the energy bands cannot 
be determined precisely, and the influence of the cluster 
size on the eigenstates of the Hamiltonian can be impor- 
tant. The differences between results obtained for a finite 
cluster, and its values for an infinite lattice, are known 
as finite size corrections. It is often difficult to extrapo- 
late the finite cluster data to the thermodynamic limit, 
and in certain cases it can lead to erroneous theoretical 
predictions. The corrections often do not decrease mono- 
tonically with the increase of the number of sites, mainly 
due to varying cluster geometries. Thus the estimation 
of finite size corrections by direct comparison of clusters 
of different sizes is difficult. Moreover, often one cannot 
compare different clusters with the same filling, as the 
number of sites as well as the number of electrons are 
integers. 

There are various methods of minimizing the finite size 
corrections. The most obvious one, the increase of the 
size of the cluster, is strongly limited by the available 
time and memory of present-day computers. However, 
there have been some recent attempts to attack this prob- 
lem, e.g., the diagonalization of the Hamiltonian in a 
reduced Hilbert space.EI Another approach to finite size 
effects is based on a specific treatment of boundary con- 
ditions. Usually, when the hopping term in the Hamil- 
tonian makes a particular jump out of the cluster, it is 
mapped back into the cluster through a translation with- 
out any change of the wave function. However, in order to 
reduce the finite size effects, twisted boundary conditions 
are sometimes used, i.e., the phase of the wave function is 
changed when the electron hops from one site to another. 
Then the properties of a larger system may be found by 
forming an average over smaller systems with different 
boundary boundary conditional In another method the 
boundary conditions are randomized by varying the mag- 
nitude rather then the phase of the "boundary" hopping. 

The aim of the present work is to evaluate the ground 
state energy of the Hubbard Hamiltonian by employing 
yet another approach to the cluster calculations. Instead 
of applying some kind of boundary condition, we propose 
to mimic the infinite lattice by treating the electron hop- 
ping into or from a cluster as a small perturbation, and 
carrying out the summation of the perturbation series. 
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II. FUNDAMENTALS OF THE 
CLUSTER—PERTURB ATIVE (CP) METHOD 

The idea of the present approach is to divide the in- 
finite lattice into small exactly soluble clusters and con- 
sider the transfer between the clusters as a perturbation 
(see Fig. @). 
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FIG. 1. The lattice divided into clusters. H c i are the clus- 
ter Hamiltonians and V is the terms describing the hopping 
between clusters. 

The Hamiltonian of the system consists of two parts: 
H = H cl + V, (1) 

where 

H d =J2&i (2) 

7=1 



is the sum of the cluster Hamiltonians, and 



(3) 



describes the hopping between nearest neighboring clus- 
ters. TV denotes the number of clusters. 

The cluster Hamiltonian Hj (1=1, ..,7V) operates on 
the Jth cluster's Hilbert space Hi (Hi : Hi — > Hi) and 
the Hilbert space of the whole system is a direct product 
of the cluster's Hilbert spaces H = Hi €3 H2 <8 •■■ €5 Hj\f. 
Taking the size of the clusters relevant for numerical di- 
agonalization, we are able to separately solve the cluster 
Hamiltonians. Now, in order to model the infinite lat- 
tice, we turn on the intercluster electronic transfer. The 
operator that moves electrons between clusters / and J 
(I and J are neighboring clusters) operates in a direct 
product of the Hilbert spaces of clusters I and J: 



Vu : Hi®Hj ^Hi®H 



(4) 



Assuming that the intercluster hopping energy is small 
compared to the distances between levels in the clus- 
ter's spectrum, we can perform perturbation calcula- 
tions, where the states of the Hamiltonian H c i will be 
the zeroth-order approximation. 



Let {|</>f)}, i = 0, 1, ... denote the set of states of 7th 
cluster, and {e^} the corresponding energy levels (eo < 
£i < £2 < ■■■)• Then 

l*{ii,i a ,..w}> = 14) ® 14) ® - ® 140 ( 5 ) 
is an eigenstatc of H c i corresponding to the energy 

E {ix,i2,-M} = e *i + e »2 + •• + ( 6 ) 

The corrections to the ground-state energy Eq in the nth 
order of the perturbation V are calculated summing the 
Goldstone diagrams 

AE n = ($ | V ( 1 v) |$ )conn, TV = 1,2, ... , 

\E n — -fi c / / 

(7) 

where |$ ) = l^o) <S> \4>o) ® ■■■ ® ) is tr± e ground state 
of the Hamiltonian H c i. 



III. APPLICATION OF THE CP FORMALISM TO 
TWO-DIMENSIONAL HUBBARD MODEL. 

On a simple square lattice, AE n vanishes for odd n 
since the perturbation V transfers an electron from one 
cluster to a neighboring one, and it needs an even number 
of jumps to return the electron to the outgoing cluster 
creating the ground state again. 

The lowest-order nonvanishing contribution to the 
ground-state energy AE2 is given by: 



A£ 2 



> \{<S>o\V\$ {ll . l2 ,.. M} )\ 2 

E - E {iui2 ,„ yi ^y 



(8) 



where the prime means summing over all states excluding 
the ground state. In the case of simple square lattice the 
symmetry of the lattice allows to simplify this expression, 
leading to 



AE 2 = 2 ^V'K^°}1^1* { hM) 
f—f 2 eo — — £i 2 



The operator V12 is given by 



Viz 



Yl (4,1. 



^04,2,0- + 0-2 l,<7 a l,2,<7 



x) +H, 



(9) 



(10) 



where we have used two indices for the creation and an- 
nihilation operators: the first indicates the position of a 
site within the cluster, and the second index is the num- 
ber of the cluster. In Fig. || the solid lines connecting 
clusters 1 and 2 represent different terms of the operator 
V12 [Eq. ©]. 
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FIG. 2. The hopping term which transfers electrons be- 
tween clusters 1 and 2 (Vu). 



The matrix element (${0.0} 
pressed in terms of (4>o\ a i <r|^i 

,1; 1 = 1, 



^12|${ii.i 2 }) Can be eX " 



and 



'o|oj )<r |^) (ct =T 
M, where Ai is the number of sites in the 



cluster), where ai^ a (a] ) creates (annihilates) an elec- 
tron with spin a on the Zth site in the cluster. Thus 
all the calculations required to evaluate the second-order 
correction are performed within Hilbert space of a size 
equal to the size of Hilbert space of a single cluster. 

In the same manner, we are able to take into account 
the fourth-order correction AE4. While the computa- 
tional potential required to perform such calculations is 
much larger (mainly due to a summation over a large 
number of intermediate states), the Hilbert space is still 
the same as in the zeroth order. The general formula for 
the fourth-order correction is too complicated to be pre- 
sented here. Instead, Fig. || shows all the diagrams that 
contribute to A £4 in the case of a simple square lattice. 
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AEf ^T A J2Y,Y1 < $ {o,o}|^2|$ { n,i 2} ) 

*i,*2 31J2 kt,k 2 
x ^{ji,j a }l^l2|^{fci,fc a })(^{fci,fc2}l^l2|^{0,0}) 

(2e - e h - e i2 )(2e - e h - e j2 )(2e - e kl - e fe ) ' 

(11) 

A contribution from each diagram is multiplied by the 
factor which reflects the symmetry of the diagram and 
the lattice; for example there are four diagrams of type 
A (obtained by rotating diagram A around site 1 by 90°), 
and each one consists of two sites, so that Ta = 2. 

Equations for contributions from other diagrams are 
more complicated, for example diagram B includes the 
following processes: 
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where I — > J describes the hopping from cluster / to 
cluster J. 

The summation over all intermediate states is a diffi- 
cult, the most time-consuming task. In order to reduce 
the computational effort we have explicitly exploited 
various symmetries of the model. For example, cal- 
culating the matrix element (^{i 1 .i 2 }\^2\^{j 1 .j 2 }), con- 
servation of the number of particles reduces the sub- 
space of states |${j lj2 j) to (Hff^-j+i O Wjv(i 2 )_i) © 
[H-Nth)-! ® Wjv(i 2 )+i) where Hn denotes a subspace of 
the cluster's Hilbert space with N electrons and N(k m ) 
is the number of electrons in the fcth state on the mth 
cluster. 
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FIG. 3. The diagrams which contribute to the fourth order 
corrections. 

The lines between clusters I and J represents a per- 
turbation Vjj that moves an electron between these clus- 
ters (in both directions, from / to J as well as from 
J to I). The interaction in Eq. ([7J) operates on the 
ground state, and as a result also produces the ground 
state. Therefore, only diagrams of the forms of loops con- 
tribute to the ground-state energy and, for example, the 
following diagram does not appear in the fourth order: 
□ —►□—»■□—>□—►□. The simplest diagram (A) gives 
the contribution which can be written explicitly as 



IV. RESULTS AND DISCUSSION 

The second- and fourth-order corrections to the 
ground-state energy were calculated for different values 
of U/t (for n = 1). The results, presented in Fig. [| and 
||, are compared with the energies obtained by exact di- 
agonalization of the Hamiltopian for 4x4 system with 
periodic boundary conditionsEj.jJ] 

Generally, apart from the region of U/t < 1, where 
the perturbation series does not converge, the results 
from both these approaches are in agreement. The cal- 
culations were performed on IBM RS/6000 workstations, 
whereas diagonalization of 4 x 4 clusters requires much 



*In Ref£3 the Hubbard Hamiltonian is written in particle- 
hole symmetric form, so the ground state energies for n — 1 
are shifted by U /4 (per site) 
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larger computing facilities (see, e.g., Ref.U). The advan- 
tage of the CP approach, comparing to the exact diag- 
onalization of larger systems is the lack of the memory 
limitations. With the increase of the order only the com- 
putational time increases, whereas the of size the diago- 
nalized matrices is constant. Of course, this method, in 
contradiction to the standard Lanczos approach, does not 
allow a calculation of the dynamical properties of a given 
Hamiltonian. The formalism can be directly applied to 
systems described by other than Hubbard Hamiltonians, 
e.g., the t — J model. 

The most attractive application of the CP formal- 
ism is the study of the hole hole effective interaction. 
Performing calculations for system with one and two 
holes in one cluster, while all the other clusters are with 
n = 1, we can calculate the two-hole binding energy: 
A = E(two holes) + E - 2E(one hole). Work in this 
direction is in progress. However, in the case of a doped 
cluster we have to perform the calculations for a degen- 
erated spectrum, where the CPU and memory usage is 
much larger. 

This scheme can be directly extended to a study of 
the nature of ground states of the undoped and doped 
systems. 



< 




2(1 




10 
U/l 



FIG. 4. The second and fourth order corrections to the 
ground state energy. 
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FIG. 5. The comparison of ground state energy calculated 
using the CP approach with the result of Lanczos diagonal- 
ization of 4 x 4 system. 
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